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Abstract 

We have studied the anisotropic three-dimensional nearest-neighbor Ising model with competitive 
interactions in an uniform longitudinal magnetic field H. The model consists of ferromagnetic 
interaction J X (J Z ) in the x(z) direction and antiferromagnetic interaction J y in the y direction. We 
have compared our calculations within a effective-field theory in clusters with four spins (EFT-4) in 
the simple cubic (sc) lattice with traditional Monte Carlo (MC) simulations. The phase diagrams in 
the h—ksT/ J x plane (h = H/ J x ) were obtained for the particular case Ai = J y / J x (^2 = Jz/ Jx) = 1 
(anisotropic sc). Our results indicate second-order frontiers for all values of H for the particular 
case A2 = (square lattice), while in case Ai = A2 = 1, we observe first- and second-order phase 
transitions in the low and high temperature limits, respectively, with presence of a tricritical point. 
Using EFT-4, a reentrant behavior at low temperature was observed in contrast with results of 
MC. 
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I. INTRODUCTION 



In recent years, the effect of a longitudinal field in the Ising antiferromagnet on an 
anisotropic simple cubic (sc) lattice has been discussed. The experimental example is com- 
pound (C 2 H 5 NH 3 ) 2 CuCl4 : The differential magnetic susceptibility x = (|f)t °^ ^ ms 
compound was analysed as a function of an extra external field (0 — 2 kOe) and of temper- 
ature (1 — 30 K; T c = 10.20 K). The compound is a typical layer-type ferromagnet, with 
a very weak antiferromagnetic coupling between the Cu 2+ layers, where has been estab- 
lished the magnetic phase diagram of the antiferromagnetic structure. One of the attractive 
points of investigating the properties of (C 2 H 5 N H 3 ) 2 CuCl^ is that as a consequence of 
the antiferromagnetic interlayer coupling we may obtain quantitative information about the 
anisotropy and Jaf (exchange coupling) by investigating the field dependence of the sus- 
ceptibility at T < T c . In previous papers [2-5] it has been reported that the Cu compounds 
of general formula {C n H 2n+ iN H^CuX^ , where n — 1, 2, 3, 4, 5, 6, 10 and X = CI or Br, 
may be considered as consisting of nearly isolated magnetic layers. Other example of the 
compound with cubic anisotropy are antiferromagnet K 2 MnF 4: jfj], (CH 3 NH 3 ) 2 MnCl4 : and 

{CD 3 ND 3 ) 2 MnCk 3- 

Three dimensional (3D) Ising models and (pseudo-) Ising physical systems have been 



analysed extensively 21^. Graim and Landau 10] studied the critical behavior of a spin- 
1/2 Ising model on a simple cubic lattice with spatially anisotropic nearest-neighbor coupling 
using the Monte Carlo method. This model on an anisotropic square lattice was investi- 
gated by using a modified mean-field theory in which the intrachain is treated exactly and 
the interactions between chains are considered in the mean- field theory llMl3| . Various ap- 
proximate methods have show this critical bahavior of the curve T/y versus H, such as mean 



field approximation (MEA) 1 151 ] . effective-field theory (EFT) [16] , mean field renormal- 
ization group (MFRG) 171]. effective-field renormalization group (EFR G) 181 ] . Monte Carlo 



simulation (MC) 19M2lj. and high-temperature series expansion (SE) 22]. For the case of 



the 3D lattice, the theoretical calculations show disagreement between differente methods. 
The results obtained by the MFA and EFT methods show a reentrant behavior in the phase 
diagram in low-temperature, i.e., if H is just above H c , then these are two phase transitions 
as the temperature is increased. 

In recent years, the effect a longitudinal field in the Ising antiferromagnetic on an 
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anisotropic square lattice was explored by MC [23]. Although MC simulations play an 
important role for the study of phase transitions and critical phenomena, the well-known 
difficulties arise when one uses standard algorithms (one-flip algorithms) [24| for the study of 



first-and second-order phase transitions. This has contributed for the deve 



opment of alter- 



native MC methods, such as parallel-tempering 25], cluster algorithms 26], multicanonical 
algorithms j^], and more recently Wang-Landau method 



In the present paper we use the MC simulations and effective-field theory in clusters 
with four spins (EFT-4). We investigate the first- and second-order phase transition in the 
plane h — k^T j J x of the Ising superantiferromagnet on an anisotropic simple cubic lattice 
in the presence of a magnetic field. Standard finite-size scaling techniques were used to 
estimate the critical temperatures. In Section II we present the model and formalism. The 
numerical results and discussions are given in Section III. Finally, the last section is devoted 
to conclusions. 



II. MODEL AND FORMALISM 
A. Hamiltonian 

The model in this work is the nearest-neighbor (nn) Ising antiferromagnetic in a longi- 
tudinal field magnetic divided into two equivalent interpenetrating sublattices A e B, that 
is described by following Hamiltonian 



n = - j x £ *?** +t + j v Y: « + £ - J. E « + t ~ H E^ C 1 ) 

■J* -I* -I* i 

l,0 X i,Oy t,O z 

where cr? is the /i(= x,y,z) component spin- 1/2 Pauli operator at site i, J x (J y ,J z ) is 
the exchange coupling along the x(y,z) axis, S x (S y ,S z ) denotes the nearest-neighbor vector 
along the x(y, z) axis and H is the longitudinal magnetic field. We define the parameters 
Ai = Jy/J x and A2 = J z /J x - The ordered state for low temperatures and fields is the 
superantiferromagnetic order (SAF), which is characterrized by a parallel spin orientation 
in horizontal direction and an antiparallel spin orientation of a parallel spin orientation of 
nearest-neighbors in vertical direction and therefore exhibit Neel order within the initial 
sub lattice A and B (see figure (CQ)). 
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If A2 = (J z = 0), the lattice is now composed of independent planes, so the model is 
exactly solved for H = 0, and the critical temperature is obtained by the following relation 
29] 



sinh ( -^k-) sinh ( ) = 1, (2) 



, k B T N J \ k B T N 

where for the particular isotropic case J x = J y = J (Ai = 1 ) we have ksT^/J = 
2/ln(l + \/2). For H 7^ 0, with Ai 7^ 1, Ai = 1 we have an Ising model with an external 



23]- 



magnetic field on an anisotropic square lattic ([I]), which was already studied by MC 
Accordingly, we improve the understanding of this model by studying it by means of the 
Effective-Field Theory and Monte Carlo simulations for the case Ai = A2 and H 7^ 0. 



B. Monte Carlo Simulation 



In order to implement the present model to perform MC by the Metropolis Algorithm, 
the simple cubic lattice of size L having L x L x L sites is decomposed into two sublattices 
(A and B) with opposite spins, corresponding to the SAF ground state. To meassure the 
SAF order, the appropriate order parameter is defined by (m s ) = ((m^ — m#)/2), where 
(m^) = /jj Ylie^ a i S ) * s the magnetization of the sublattice fi = A, B and N = L 3 number 
of spins. The susceptibility related to this order parameter is defined as follows: 



X = L d ((ml) - (m s ) 2 )/T, 



(3) 



where d is the dimension of the lattice. 

In our simulations we have considered lattices with periodic boundary conditions. In 
order to determine the system's behavior in the thermodynamic limit (L — > 00), which is 
imposible to implement on account of computational limitation, we have to use a finite-size 
scaling theory 30[. Accordingly, to locate the critical temperature for second-order phase 
transitions, we approximately locate the crossing point of curves for different sizes of the 
fourth-order cumulant U^(L) (Binder Cumulant) defined as [31 ] 



U A {L) 



3K> : 



(4) 
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where (m 2 s ) and (Trig) are the canonical averages of the second and fourth moments of mag- 
netization, respectively. On the other hand, the critical temperature T c can also be obtained 
by means of the relation T C L = T c + aL^ 1 ^, where T C L is the pseudocritical temperature 
corresponding to the size L, and v is the correlation critical exponent. For first-order phase 
transitions the same formula applies by setting \ jv = d. A more extensive description of 



the Monte Carlo method was published elsewhere |20j and the reader is referred there for 
further details. We performed simulations for Xi = X 2 , for several values of h, with a tem- 
perature step AT = 0.01 and runs comprising up to 2 x 10 5 MC steps after equilibration. 
The statistical errors of the MC simulations used for the estimation of T^(Xi, X 2 ,h) of a 
particular Ai, A2 and h = H/J x were found much smaller than the statistical errors coming 
from the fact that we used. Therefore, the error bars are not shown in our graphs because 
they are smaller than the symbol sizes. 

C. Effective-Field Theory 

As a starting point, the averages of a general function involving spin operator components 



0({n}) are obtained by [32] 



(0{{n})) = { Tr{n}e -,n {n} ) . (5) 

where the partial trace Tr{ ra } is taken over the set {n} of spin variables (finite cluster) 
specified by the multisite spin Hamiltonian H{ n y and (• • ■) indicates the usual canonical 
thermal average. 

The method treats the effects of the surrounding spins of a finite cluster with N spins 



through a convenient differential operator technique |33] such that, in contrast to the usual 
MFA procedure, all relevant self-spin correlations are taken exactly into account. The inter- 
actions within the cluster are exactly treated and the effect of the remaining lattice spins is 
treated by a given approximation (here we use the random phase approximation- RPA). 

To treat the model (JT|) by the EFT approach, we consider a simple example in cluster of 
size N = 4 spins, and the Hamiltonian for this cluster is given by 

4 

- (3Hi = -KX x o{a z 2 + KX 2 o z 2 a\ - KX^at + KX 2 a\a{ + ^ a r a z r , (6) 

r=l 
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where a r = L — KJ2s r a r+Sf, with K = (3J X , Ai j2 = J y , z / Jxi L = (3H and 5 r corresponds to 
nn vectors. 

Substituting Eq. (EJ) in (j5J), we obtain the average magnetizations in sublattices A and 
B, respectively, by 

/ din 24(a) \ 

mA = fa) = ( w7w-\ > , (7) 



and 



, / 9 In 24(a) . 



with 

Z 4 (a) = Tr w e-^ 4 , (9) 

where a = (ai, a 2 , a 3 , a 4 ) and {a} = {er*, erf, erf, cr|}. 

Using the identity exp(a • D)f(x) = f(x + a), where D = (Dx, D 2 , D 3 , D 4 ) and x = 
(xi, X2, Xs, x 4 ) are four- dimensional differential operator and vector, respectively, = 
and also the van der Waerden identity for a* component Pauli spin operator, i.e., exp(Aof) = 
cosh(A) + o\ sinh(A), Eqs. (J7J) and (jSJ) are rewritten as [p, — A or B) 



with 



and 



lz-2 z-2 

m » = { n («i + °ia+sA) n («a + °2b+&.a) 

Ul ,- a 52 v (10) 

II («3 + OaB+is^a) 11 («4 + ^+54^) ) /m( x + L )lx=0 

<91nZ 4 (x) ^1 + ^2 + ^3 + ^4 

<9xi 01 + 02 + 03 + 04 

<91nZ 4 (x) ^1-^2 + ^3-^4 , 1Q v 
/M X J - — 7 — — —7 — —7 — — — ) { LZ 

OX 2 01 + 02 + 03 + 04 

where ip x = sinh(C 1 ) + e 2K[ - Xl+x ^ sinh (C 2 ), ^2 = sinh (C 3 ) + e 2X ( Al ~ A2 ) sinh (C 4 ), 
^ 3 = sinh (C 5 ) + e 2X (- Al+A2 ) sinh (C 6 ), ^4 = sinh (C 7 ) + e - 2 ^ Al+A2 ) sinh (C 8 ), X = 
cosh(Ci) + e 2/f ( Al + A2 )cosh(C 2 ), 2 = cosh(C 3 ) + e 2X ( Al ~ A2 ) cosh (C 4 ), 3 = cosh (C 5 ) + 
e 2K( " Al+A2) cosh(C 6 ), 4 = cosh(C 7 ) + e" 2 ^ 1 ^ 2 ) cosh (C 8 ), C x = x x + x 2 + x 3 + x 4 , 
C 2 = xi -x 2 -x 3 + x 4 , C 3 = 2i+x 2 -x 3 + x 4 , C 4 = xi -x 2 + x 3 + x 4 , C 5 = x x -x 2 + x 3 —x 4 , 
Cq = x\ + x 2 + x 3 — £4, C7 = x 4 — X2 — x 3 — X4 and C$ = x x + x 2 — x 3 — X4. 

The magnetization in Eq. (jTUl) is expressed in terms of multiple spin correlation 
functions. The problem becomes unmanageable when we try to treat exactly all boundary 



spin-spin correlation function present in Eq. (fit)]) . Here we use a decoupling procedure that 
ignores all higher-order spin correlations on both right-hand sides in Eq. f JTOj) . namely 



'^iA^jB ■ ■ ■ °\Al - m Am B ...m A , 



'lAl — ni'Aili'B ■ ■ ■ "<V1; (13) 

where i ^ j ^ ■■■ ^ I and m M = (c? ) (/i = A, B). The approximation ( fl3l) neglects 
correlations between different spins but takes relations such as ((of,,) 2 ) = 1 exactly into 
account, while in the usual MFA all the self- and multi spin correlations are neglected. We 
can then rewrite the Eq. f JT0|) in the form 



m A = ("i + m A Pif (a 2 - m B (3 2 ) 2 (a 3 - m B (3 3 ) 2 (a 4 + m A (3 4 ) 2 f A (x) \ , (14) 



and the expression for the magnetization in sublattice B is given by 



m B = (at + m B f3 1 ) 2 (a 2 - m A (3 2 ) 2 (a 3 - m A f3 3 ) 2 (a 4 + m B (3 4 ) 2 /b(x)| 



x=o • 



(15) 



Defining the uniform m = ^(m A + m B ) and staggered m s = \{m A — m B ) magnetizations, 
and using the identity exp (a • D) F(x)| x = -F(a), we obtain 



m s = A(m s , m, T, H) = A 2r+ i(m, T, H)m 



2r+l 



(16) 



r=0 



and 



m 



J2B 2r (m,T,H) 



2r 



(17) 



r=0 



where the expressions for the coefficients A p (m,T, H) and B p (m,T, H) are again omitted 
here. 

We note that is not possible to calculate the firts-order transition line on the basis of 
only the equation of state, Eq. ( 1T6|) . to solve this problem one needs to calculate the free 
energy for each state (P and SAF). Assuming that this equantion of state is obtained by 
the minimization of a given free energy functional like $(m s ) (i.e., 5$ = 0), then after 
intergration we obtain 



$(m s ) = A 1 (T,#) + A 2 (T,#) 



-X)A 2r+1 (m,T, J ff) 



2r+l 



r=0 



2r + 2 



(18) 



S 



where Ai^T, H) are arbitrary functions which turn out to be irrelevant for searching the 
second and first-order transitions. 

To obtain this phase transitions we use Maxwell construction, that correspond to the 
intersection point where the free energy between the phases are equal. In the case of the 
transitions between the SAF ordered and P disordered (m s ) phases we obtain the point of 
intersection $sAF(ms) = $p(m s ) from Eq. (fl8|) . i.e, 



The phase transition temperatures between the P and SAF states are found by simultane- 
ously solving three transcendental expressions, Eqs. ( |T6i) . ( IT71) and ( |T9i) . For a second-order 
transition, we obtain m s = 0, while the first-order transition we have m s 7^ 0, where this 
value corresponds to the discontinuity of the staggered magnetization at T = T*(H). 

III. RESULTS AND DISCUSSION 

We obtained the phase diagrams in the h — ksT j J x plane for Ai = A2 = 1 of the model [T] 
by using EFT-4 and MC simulatons. Both methods confirmed the existence of a tricritical 
point in a critical frontier separating the SAF order with the paramagnetic disorder as 
shown in Figure 2. For completeness, Figure 2 also includes the critical frontier obtained 
by MC simulatons for the particular case Ai = 1, A2 = (square lattice), which lacks of 
first-order criticality j^. In what effective-field calculations concerns, we obtain the critical 
frontier, which consists of the transition temperature as a function of the external magnetic 
field, separating the SAF and P phases, by simultaneously solving the three transcendental 
equations Eqs. ( 1T61) . ( |T71) and ( TT8i) . The range of ratio h = H/J x determines the limits of 
second-order (0 < h < h t ) and first-order (h t < h ^ 2) phase-transition frontiers, where 
(h t , kg/ J x T t ) is the tricritical point (TCP). From Landau theory, the TCP is calculated 
by the condition given by Ai(m,T t , H t ) = 1, A 3 (m,T t , H t ) = and A 5 (m,T t , H t ) > 0, 
obtaining in EFT-4 the following values: h t = 1.72 and ksTt/Jx = 3.42. The temperature 
variation of m s (T,H) for fixed reduced field H/J x present two type of behavior. The first 
for low-field the order parameter decrease to zero continuously as the reduced temperature 
approaches to the critical point, the temperature (H fixed) where m s = is the second-order 
phase transition temperatute, T C (H). On the other hand, when we are at high-fied the 




r=0 



(19) 
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temperature at which the order parameter make a jump discontinuity is the first-order 
phase transitions going to zero discontinuously at the point T C (H). For h value greater than 
the upper limit of these field (i. e., h > h c = 2.0), the system exhibits no phase transition 
(the order parameter is null for all finite temperature, T > 0) 



In order to obtain relevant critical points to get the critical frontier in MC simulations, 
we estimate the critical temperature for a given value of field by a finite-size scaling 
analysis. For instance, for h = H/ J x = 1.1, Figure 3 shows how the critical temperature was 
located around the crossing point of U± — k B T/J x curves for sizes L = 16,32,64, resulting 
in k^TcjJx ~ 4.14. Note that critical temperatures for lower fields are superestimated 
by the EFT-4 technique (see Figure 2). The limit of the second-order frontier was found 
by obtaining several U4 — kBT/J x curves for different field values for a given lattice size. 
Consequently, we found the beginning of the first-order criticality after detecting a change 
in the behavior of the Binder Cumulant curve above a certain value of field. This is well 
exemplified in Figure 4, for two different field values h = 1.84, 1.88, for L = 30. Figure 
4a shows that the Binder Cumulant curve already presents a clear first-order transition 
for h = 1.88, because of the sharply negative fluctuated peaks around the pseudocritical 
temperature. Moreover, Figure 4b also suggests this first-order criticality by the aparent 
Delta-Dirac form of the susceptivility, for h — 1.88. In order to confirm what is suggested 
in short lattice sizes, we draw upon the fact that around the critical temperature of a 
phase transition the susceptibility peaks behave as Xmax ~ L e , where 9 — d (d is the 
lattice dimension), for a first-order criticality, whereas 9 — j/u for a second-order one. 
Consequently, Figure 5 shows a finite-size scaling of the susceptibility peak versus L, for 
sizes L = 15, 20, 30 and 40. For h = 1.84, we estimated 9 = 1.98(3), which is (within the 



error bar) an universal value 



34 



35] for the second-order criticality of the three-dimensional 
Ising model. For h = 1.88, 9 ~ 2.8, which is close to d — 3, which suggests a first-order 
transition. Therefore, the tricritical point must be in the interval 1.84 < h < 1.88. On the 
other hand, the first-order curve must end at h c — 2. 

We could estimate some points of the first-order frontier by a finite-size scaling method, 
however, in order to avoid finite-size problems to obtain well behavied curves for the specific 
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heat or the susceptibility around critical points belonging to the first-order frontier, we 
decided to estimate the corresponding critical temperatures for given fields by computing 
the free energy versus (3 = l//cgT, from Monte Carlo data of energy versus /3, for a given 
size. For instance, in figure 6a it is shown the hysteresis effect from energy data around 
the critical temperature obtained when cooling or heating the system too fast. This is 
because of huge intrinsic autocorrelation times in finite-size systems. We can improve the 
location of the critical temperature by obtaining the associated free energies of the low- 



and high-temperature branches 



361 ] . Consequently, the crossing point of the free energy 



branches gives a good estimation of the critical temperature, inasmuch the stable phase has 
the lower free energy, as shown in figure 4b. 



The first-order frontier in MC simulations does not present a reentrant behavior as 
EFT-4 frontier does. So, this reentrance seems to be an artifact of the effective-field 
approach. 



IV. CONCLUSIONS 



In summary, we investigated the anisotropic three-dimensional nearest-neighbor Ising 
model with competitive interations in an uniform longitudinal field. We obtained the phase 
diagram in the plane h — ksT/J x , where the critical frontier separates the SAF order with 
the paramagnetic disorder. For lower fields, the EFT-4 approach superestimates the critical 
temperatures obtained by MC simulations. A tricritical point was found for 1.84 < h < 1.88, 
by a finite-size analysis. At low temperatures EFT-4 calculations show a reentrant first- 
order frontier, which does not appear by MC simulations. It suggests that improvements in 
treating correlations, or by increasing the cluster size in the Effective-Field approach, could 
correct this reentrant curve. At zero temperature, the critical field is exactly obtained, so 
h c = 2.0. Our quantitative estimation for the tricritical point could be bettered by using 
larger lattice sizes with better Metropolis Monte Carlo techniques like Parallel Tempering. 
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FIG. 1: Configuration of the super antiferromagnetic system of spins used to treat the model whose 
hamiltonian is described in Eq. ([!]). 
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FIG. 2: Critical frontiers in the plane h — ksT/Jx (h = H/ J x ). These frontiers separate the colinear 
order (SAF) and the paramagnetic phase (P). The curve (a) was obtianed by MC simulations 
231 ] for the present model implemented in square lattices. The solid line corresponds to second- 
order transitions. The curves (b) and (c) are our results obtainded by MC and EFT-4 methods, 
respectively. The dashed lines corresponds to first-order transitions appeared for the present model 
implemented in cubic lattices. 
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FIG. 3: The inset shows the fourth-order cumulant U^L) for a particular field h = H/J x = 1.1 
and system sizes L = 16, 32 and 64. 
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FIG. 4: The two figures show, two different phase transition orders at different values of the 
external field, through the Binder Cumulant (a), and through the susceptibility of the relevant 
order parameter (b), for the present model implemented in the sc for L = 30. 
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FIG. 5: Finite-size scaling of the susceptibility related to the relevant order parameter for two 
values of field, corresponding to the present model implemented in the cubic lattice. The different 
estimated slopes suggest a tricritical point for 1.84 < h < 1.88. 
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FIG. 6: (a) Heating and cooling Monte Carlo runs in (3 = 1/ksT, for the energy, (b) The associated 
free energy obtianed from the Monte Carlo data in (a), by thermodynamic integration. 
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